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1. Introduction 

Hydrodynamical models have certain advantages over transport model 
calculations in describing heavy ion collisions. One of the most important 
is that, once the equation of state and initial conditions of the matter are 
specified, the evolution of the system is determined. No knowledge of the 
underlying microscopic processes is required. This is especially important 
when studying the predicted phase transition from hadronic to partonic 
^ ■ degrees of freedom (and vice versa) - a process for which details are still 

unknown. 

The hydrodynamical description is relatively simple and fulfills the con- 
servation laws without additional constraints. The use of familiar concepts 
like temperature, pressure and flow velocity also provides an intuitive and 
transparent picture of the evolution. The price to be paid for these advan- 
tages is a set of bold assumptions: local kinetic and chemical equilibrium 
and lack of dissipation. This set of assumptions may or may not be valid in 
such a small system as that formed in a heavy ion collision. 

In a hydrodynamical description the evolution is assumed to proceed as 
follows: In the initial collision a large fraction of the kinetic energy of the 
colliding nuclei is used to create many secondary particles in a small volume. 
These particles will collide with each other sufficiently often to reach a state 
of local thermal equilibrium. When the system has reached local equilibrium 
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it is characterized by the fields of temperature, T(x), chemical potentials 
associated with conserved charges, Hi(x), and flow velocity, u^(x). The 
evolution of these fields is then determined by the hydrodynamical equations 
of motion until the system is so dilute that the assumption of local thermal 
equilibrium breaks down and the particles begin to behave as free particles 
instead. 

The numerical solution of hydrodynamical equations of motion in all 
three spatial dimensions is a tedious problem. In most approaches some ap- 
proximate symmetry is applied to reduce the number of spatial dimensions 
where numerical solution is needed to two or one. 

In the so-called Bjorken model [1] the main idea is the boost invariance of 
the longitudinal flow. The longitudinal flow is assumed to be given by v z = 
z/t at all times. This leads to particularly simple solutions of the equations 
of motion since the longitudinal expansion can be solved analytically. Also 
it is sufficient to solve the equations of motion in the transverse plane at 
z = since the solution is independent of the boosts along the beam axis. 
The obvious drawback in this approximation is that the observables are 
independent of rapidity. 

In central collisions of spherical nuclei the expansion can be simplified 
by assuming cylindrical symmetry. Of course this symmetry can not be 
applied to non-central collisions where the shape of the source has a crucial 
role in the buildup of elliptic anisotropy. 

2. Hydrodynamical models 

2.1. The basics 1 

Hydrodynamics is basically an application of conservation laws. The 
local conservation of energy and momentum and any other conserved four- 
currents , i = 1, . . . , n are expressed by 

dfjT^ = and drf? = 0, 

respectively, where T^ u is energy momentum tensor. Without any addi- 
tional constraints these 4 + n (n is the number of conserved currents) equa- 
tions contain 10 + 4n unknown variables. The simplest and most commonly 
used approach to close this system of equations is the ideal fluid approxi- 
mation which reduces the number of unknown variables to 5 + n. 

In the ideal fluid approximation the energy momentum tensor of the 
kinetic theory, 

1 For a more detailed discussion see ref. [2]. 
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and currents ji are supposed to have forms 

T^ u = (e + v)^u v - pg^ and jf = mvP, 

where e, p and rii are energy density, pressure and number density of charge 
i in the local rest frame of the fluid, and is the flow four-velocity of 
the fluid. In other words all dissipative effects, such as viscosity and heat 
conductivity, are assumed to be zero and the fluid is always in perfect local 
kinetic equilibrium. The additional equation needed to close the system 
of equations is provided by the equilibrium equation of state (EoS) of the 
matter, which connects the pressure to the densities: P = P(e,rii). 

In principle it is possible to include small deviations from local ther- 
mal equilibrium by including dissipative effects, but in practice relativistic 
viscous hydrodynamics is very difficult to implement and has not yet been 
done [2]. For preliminary results and estimates of the effects of viscosity, 
see ref. [3] . 



2.2. Equation of state 

The present results from lattice QCD calculations point to a phase 
transition from hadronic to partonic degrees of freedom at a temperature 
T c ~ 155 - 175 MeV, but the order of the phase transition is still uncer- 
tain [4]. So far no nuclear equation of state (EoS) based on lattice results 
has been employed in hydrodynamical calculations, mostly because lattice 
calculations are available only at zero net baryon density. 

The usual way to construct an EoS is to use the EoS of a hadron gas at 
low temperatures and the EoS of an ideal parton gas with a bag constant 
at temperatures above the critical temperature T c (see e.g. ref. [5]). The 
EoS of an interacting hadron gas is approximated by the EoS of an ideal 
resonance gas with resonances up to 1.5 - 2 GeV mass. It is known that the 
inclusion of higher- lying resonances mimics interactions between hadrons 
well in temperatures up to the pion mass [6] , but there is no reliable way to 
check whether this holds at higher temperatures [5]. The phase boundary 
is determined by using the Gibbs criteria. A first order phase transition 
between the hadronic and partonic phases is achieved by connecting the 
EoSs with the Maxwell construction. 

This procedure is thermodynamically consistent and the EoS both above 
and below the phase transition temperature is based on well established 
models. However, one of its disadvantages is that it is not possible to 
determine the phase transition temperature and latent heat independently. 
To circumvent this drawback Teaney et al. [7] took only the speed of sound, 
c s = i/dp/de, from the bag model EoS and made the critical temperature 
and the latent heat explicit parameters of their model. 
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Constructing an EoS with a second order phase transition or crossover 
between the phases when only the EoSs of the separate phases are known 
is nontrivial. To achieve this in a consistent way Zschiesche et al. [8] con- 
structed a family of EoSs based on a parametrized a — to model. By chang- 
ing the values of the parameters they were able to create two EoSs with 
a first order phase transition but different latent heats and an EoS with a 
crossover phase transition. Strictly speaking these EoSs do not contain a 
deconfinement phase transition but a chiral phase transition. However the 
EoS below and above phase transition temperature is very similar to the 
more conventional constructions explained above. 

2.3. Initialization 

Local thermal equilibrium is one of the assumptions of a hydrodynami- 
cal model; the model itself does not specify the mechanism that leads to an 
equilibrated state. Since at RHIC energies the initial particle production 
is definitely not an adiabatic process, hydrodynamics can not be used to 
describe the initial collision. The hydro dynamical evolution must begin at 
a sufficient time after the initial collision when the system has had time 
to reach thermal equilibrium. The initial state of the system, i.e. the den- 
sity distributions and flow velocities at the beginning of the hydrodynamic 
evolution, are not given by the model either but must be given as external 
input. 

When a boost invariant expansion is assumed, the choice of an initial 
state is reduced to a choice of transverse density and velocity profiles. A 
simple approach is to fix the value of the entropy to reproduce the observed 
final particle multiplicity and to distribute it on the transverse plane as- 
suming a constant density profile within the radius of the colliding nuclei 
using a Fermi function to smooth the edges of the system. However, this 
approach cannot be applied to non-central collisions. One must also remem- 
ber that since it is the local pressure gradients which drive the development 
of transverse flow and the evolution of the system, the details of flow are 
sensitive to the details of the initial distributions [9]. In the same way 
the final anisotropics are proportional to the deformation of the source and 
additional constraints to the initial distributions are required. 

It is known that up to SPS energies the multiplicity scales with the num- 
ber of nucleons participating in the collision [10]. On the other hand, in the 
high energy limit one expects the individual parton-parton collisions to con- 
tribute equally to primary particle production and therefore the multiplicity 
should scale with the number of binary collisions [11]. Therefore it is nat- 
ural to initialize the system using a localized version of these approaches: 
to assume that the density is proportional either to the number of partici- 
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pants or to the number of binary collisions per unit area in the transverse 
plane. Both approaches, or a combination of them can be used to fix the 
initial entropy or energy density; a comparison was made in ref. [12]. The 
initial transverse flow velocity is customarily assumed to be zero, although 
pre-equilibrium density gradients might lead to small, but finite, transverse 
flow velocities at the time of thermalization. 

If the assumption of boost invariance is relaxed the choice of initial 
state becomes considerably more complicated. There are few constraints 
for the flow velocity profile or the longitudinal density distributions. The 
choice of a particular parametrization and the values of the parameters is 
largely based on trial and error - tuning the model until a reasonable fit to 
experimental rapidity distributions is achieved. Even for the same EoS there 
are several possible initial states which lead to an acceptable reproduction 
of the data [16]. For a sample of initial profiles used successfully see refs. [5, 
13, 14, 15]. 

An alternative approach to determine the initial state is to use some 
other model to calculate it. For example, event generators [17] or perturba- 
tive QCD (pQCD) calculations [18] have been used for this purpose. Even 
if these approaches increase the predictive power of hydrodynamics, ther- 
malization is still an additional assumption. 

2.4- Freeze- out 

At some point in the evolution particles will begin to behave as free 
particles instead of a fluid and the hydrodynamical description breaks down. 
When and where that happens is not given by hydrodynamics but must be 
included as an external input. The conventional approach is to assume this 
to take place as a sudden transition from local thermal equilibrium to free 
streaming when the mean free path of the particles becomes larger than the 
system size, or when the expansion rate of the system is larger than the 
collision rate between particles. Finding where these conditions are fulfilled 
is a nontrivial problem. Since the mean free path is strongly dependent 
on temperature the usual approximation assumes that the freeze-out takes 
place on a hypersurface where temperature (or energy density) has a chosen 
freeze-out value. This temperature is of the order of the pion mass, but its 
exact value is largely a free parameter which can be chosen to fit the data. 
In Pb+Pb collisions at the SPS the values of freeze-out temperatures vary 
between 100 and 140 MeV in different calculations [19]. 

These values are somewhat smaller than the T ~ 160 MeV freeze-out 
temperatures obtained using thermal models to fit the particle abundan- 
cies at final state [20]. This can be understood by noticing that the former 
approaches assume kinetic equilibrium while the thermal models assume 
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chemical equilibrium. Since chemical equilibrium requires frequent inelastic 
collisions, while kinetic equilibrium only requires elastic collisions, it is nat- 
ural to assume that inelastic collisions cease first and chemical freeze-out 
occurs at higher temperatures than kinetic freeze-out. Thus the system may 
be in local kinetic, but not chemical, equilibrium at the later stages of its 
evolution. How much this would affect the EoS and whether this change in 
the EoS would have any observable effects in the evolution of the system is so 
far largely unexplored (with some early exceptions like ref. [21]). However, 
first preliminary results have been shown and more are in preparation [22] . 

After choosing the surface where the freeze-out takes place, the thermo- 
dynamic variables characterizing the state of the fluid must be converted to 
spectra of observable particles. A practical way of doing this is the Cooper- 
Frye algorithm [23] where the invariant momentum distribution of a hadron 
h is given by 



Here the temperature T(x), chemical potential fj,(x) and flow velocity u^{x) 
are the values on the decoupling surface ct/. Besides its relative simplicity, 
this approach has the advantage that if the same equation of state is used on 
both sides of decoupling surface, both energy and momentum are conserved. 
However, the Cooper-Frye formula has a conceptual problem. At those areas 
where the freeze-out surface is spacelike, the product p^da^ may be either 
positive or negative, depending on the value and direction of p^. In other 
words, the number of particles freezing out on some parts of the freeze-out 
surface may be negative. These negative contributions are small (few per 
cent, see ref. [7]) and usually ignored. More refined procedures without 
negative contributions have been suggested [24] but their implementation 
is complicated. A model using one of these refined freeze-out procedures 
is in preparation [25] and it will be interesting to see how large an effect 
the freeze-out procedure has on the final particle spectra in a full-fledged 
calculation. 

Another way to refine the hydrodynamical freeze-out procedure is to 
switch from a hydrodynamical to a microscopic transport model description 
well within the region where hydrodynamics is supposed to be applicable [7, 
26]. Besides giving a better description of freeze-out, such models include 
the separate chemical and kinetic freeze-outs. The main drawback of such 
models - apart from the increased complexity - is that the region where the 
switch from hydro to transport description should take place is as poorly 
defined as kinetic freeze-out surface in ordinary hydrodynamical calculation. 
The educated guess employed in both refs. [7, 26] is that the switch happens 
immediately after hadronization. 
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Fig. 1. Charged particle multiplicity 
calculated at different collision energies 
and compared to Phobos data [27] and 
calculations (solid and dotted lines) in 
ref. [29]. A=197 and 208 refer to the 
mass of the nuclei and A e s to a finite 
impact parameter. Figure and calcula- 
tion are from ref. [18]. 
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Fig. 2. Charged pion and antiproton 
spectra from most central Au+Au col- 
lisions at the y/sNN — 130 GeV as mea- 
sured by the Phenix [31] and STAR [32] 
collaborations. The hydrodynamical 
calculation and figure are from ref. [33] . 



3. Comparison with the data 

3.1. Charged particle multiplicity and transverse energy 

In ideal fluid hydrodynamics the final particle multiplicity is propor- 
tional to the entropy of the initial state. Thus in a certain sense multiplicity 
is not really a hydrodynamical result, but something one dials in as input. 
Nevertheless, when the initial state is not chosen by comparing the final 
result to data, but calculated from perturbative QCD as in ref. [18], the 
final particle multiplicity becomes a prediction of the model. 

In ref. [18] Eskola et al. used the EKRT saturation model to calculate 
the production of minijets in the primary collisions and then converted 
this result to the initial state of hydrodynamic evolution. Their results 
for charged particle multiplicities at RHIC and LHC energies are shown 
in fig. 1. The results agree well with the data measured by the Phobos 
collaboration [27, 28]. It is important to remember that no fitting or fine 
tuning was done to achieve this result, and that the multiplicity at ^/snn = 
200 GeV energy was not yet measured when this calculation was done. 

When compared with the EKRT saturation results shown in ref. [29], 
which did not contain an expansion stage, one notices that hydrodynamical 
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expansion causes only a small change in multiplicity, but the transverse 
energy decreases by a factor 3. However, at RHIC the calculated E t is still 
10-20 % larger than the experimental value [30]. 

3.2. pt spectra 

The conventional way to initialize a hydrodynamical calculation is to 
use experimental hadron data to fix the initial values. Depending on the 
details of the initialization the calculated pt distributions can be only fits 
to the data or have some predictive power. In ref. [33] the initialization is 
done by first choosing the initial entropy and the net baryon densities to 
reproduce the observed pion multiplicity and p/p ratio in the most central 
collisions. Then a combination of different parametrizations mentioned in 
section 2.3 and ref. [12] is chosen to reproduce the observed multiplicity 
per participant as function of centrality. In this process the slopes of the 
spectra or the centrality dependence of the slopes are not used as an input 
and their calculated values can be taken as predictions. 

In central collisions this initialization led to a maximum initial temper- 
ature T max = 328 MeV and an energy density e max = 21.4 GeV/fm 3 at an 
initial time To = 0.6 fm/c. At time r = 1 fm/c, usually used to estimate 
initial energy density via Bjorken's formula [1], the corresponding average 
energy density is (e) = 5.4 GeV/fm 3 , which is consistent with experimental 
estimates [30]. 

Local chemical equilibrium is assumed to hold until kinetic freeze-out. 
Thus chemical and kinetic freeze-outs take place at the same temperature. 
Even if the pion yield is correctly reproduced in ref. [33], it is not possible 
to obtain both correct proton and antiproton yields simultaneously. The 
authors have chosen to circumvent this problem by calculating the particle 
yields at hadronization temperature T c = T c h em = 165 MeV and calculating 
the slopes of the pt spectra at Tf = 128 MeV. Subsequently they rescale the 
particle yields at kinetic freeze-out to their values at chemical freeze-out by 
hand. The results shown here in figs. 2 and 3 are very similar to those in 
ref. [7]. Since separate chemical and kinetic freeze-out was included in the 
transport model description of hadronic phase of ref. [7], one can conclude 
that the method applied in ref. [33] provides a reasonable approximation. 

It has been observed that around pt = 2 GeV the antiproton and pion 
yields are roughly equal, and the slopes of the distributions suggest that 
the antiproton yield is larger than pion yield for p t > 2 GeV [31]. As can 
be seen in fig. 2 this phenomenon can be explained as a simple consequence 
of strong transverse flow. Besides giving an explanation for this so-called 
anomalous p/ir ratio, hydrodynamical calculation provides very good fit to 
both pion and antiproton spectra both in central and semicentral collisions. 
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Fig. 3. Charged pion (left panel) and antiproton (right panel) spectra from Au+Au 



collisions at ^snn = 130 GeV at various centralities as measured by the Phenix [31] 
and STAR [32] collaborations. The hydrodynamical calculation and figures are 
from ref. [33]. 

As can be seen in fig. 3 the deviation from the data is significant only in 
the most peripheral collisions. 

3.3. Elliptic anisotropy 

Since the initial particle production is azimuthally symmetric, azimuthal 
anisotropy of the final particle distributions is a signal of rescatterings 
among produced particles. More frequent rescattering can be expected to 
lead to a larger anisotropy and since hydrodynamics assumes zero mean free 
path and thus an infinite scattering rate, it provides an upper limit to ob- 
servable anisotropics. Anisotropy is quantified by measuring the harmonic 
coefficients v n (y,pt; b) of a Fourier expansion in <p p of the measured hadron 
spectrum dN / (dy pt dpt d<p p ) [34]. Anisotropy characterized by a non-zero 
second coefficient, v%, is called elliptic anisotropy or elliptic flow [35]. 

In fig. 4 the centrality dependence of the elliptic anisotropy coefficient 
i>2 in Au+Au collisions at ^snn = 130 GeV energy is shown [36]. The 
calculations of refs. [7, 37] give very similar results: in central and semi- 
central collisions the data reaches the hydrodynamical limit and depending 
on the EoS and freeze-out temperature (the latter is not shown in fig. 4) 
the agreement is satisfactory even close to peripheral collisions. The cal- 
culation shown in fig. 4 was carried out using various EoSs with different 
latent heats [7]. The best fit to the data is obtained using an EoS with 
a relatively large latent heat, but when other observables are considered, 
the authors conclude that LH8 EoS with smaller latent heat leads to best 
overall description of data. 

In a hydrodynamical model the final anisotropy is proportional to the 
initial deformation of the source. The deformation as a function of impact 
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Fig. 4. Elliptic anisotropy coefficient V2 for charged particles at RHIC as a function 
of the number of participants (relative to the maximum) as measured by the STAR 
collaboration [36] and calculated using different equations of state [7], The label 
RG corresponds to a resonance gas EoS without a phase transition. The number 
in other labels describes the latent heat. The figure is from ref. [7], 



parameter depends on the particular parametrization applied but, as shown 
in ref. [12], when the anisotropy is averaged over all momenta and central- 
ity is expressed as a fraction of the total multiplicity, differences between 
parametrizations even out. Thus the result shown in fig. 4 is independent 
of the particular parametrization. 

Thept differential anisotropy, V2(pt), for pions and for the sum of protons 
and antiprotons [38] depicted in fig. 5 also shows clear hydrodynamical 
behaviour. As predicted 2 in ref. [39], the heavier the particle, the smaller 
the anisotropy at low p t . So far there are no published data regarding the pt 
differential anisotropy of strange particles so whether they also obey this rule 
remains to be seen. The identified particle data is so far limited to the low pt 
region shown in fig. 5. The pt differential anisotropy of negative hadrons [40] 
follows the hydrodynamical calculations up to transverse momenta pt < 1.5 
- 2 GeV (not shown) where the anisotropy saturates. The deviation can be 
understood as a sign of incomplete thermalization of high-pt particles. 

It is well known that the changes in the EoS can be compensated by 
changing the initial state and freeze-out temperature (see e.g. ref. [16]). A 
remarkable feature of fig. 5 is that changes in the EoS affect the anisotropy 
of both pions and nucleons in the same way - a stiffer EoS leads to a larger 

2 A similar result was later obtained in ref. [7]. 
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Fig. 5. p t differential elliptic anisotropy of pions (left panel) and pro- 
tons+antiprotons (right panel) in minimum bias collisions as measured by the 
STAR collaboration [38] and calculated using different equations of state and freeze- 
out temperatures [39]. The letters Q and H in the labels stand for an EoS with a 
first order phase transition and a hadron gas EoS without a phase transition, re- 
spectively. Numbers in parentheses stand for the freeze-out temperature in MeV. 

anisotropy at low values of pt - but a change in freeze-out temperature 
changes the pion and nucleon anisotropies in opposite directions. A lower 
freeze-out temperature leads to a larger anisotropy for pions but to a smaller 
anisotropy for nucleons (for further discussion, see ref. [39]). This may 
provide an additional method of constraining possible equations of state 
and freeze-out temperatures. 

In the studies mentioned above the expansion was assumed to be boost 
invariant. Judging by the measured rapidity distributions [41], this is a rea- 
sonable assumption close to midrapidity, but it makes it impossible to make 
any statements about the rapidity dependence of any variable. To study 
the rapidity dependence of the elliptic anisotropy the assumption of boost 
invariance was relaxed and the calculation was done using a genuinely three 
dimensional model in ref. [15]. When compared to the excellent agreement 
with the data in figs. 4 and 5, the result depicted in fig. 6 [15] may look less 
satisfactory. The data [42] reaches the hydro dynamical value only around 
midrapidity. On the other hand, even this result reproduces the data within 
a one to two units of pseudorapidity wide window. This area already con- 
tains most of the produced particles. It is also worth remembering that 
anisotropy in hydrodynamical models depends strongly on the initial shape 
of the system. The initialization used in ref. [15] is relatively simple and 
more sophisticated initialization may lead to better fit to the data. There- 
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Fig. 6. Pseudorapidity dependence of elliptic anisotropy for charged particles in 
minimum bias collisions at -^saT/v = 130 GeV energy. The figure and calculations 
are from ref. [15] and the experimental data from ref. [42]. Dashed and dotted lines 
correspond to different initializations of the model, see ref. [15]. 

fore it is premature to conclude based on this data and calculation alone 
that thermalization is reached only at midrapidity 

3.4- Two particle correlations 

Two particle momenta correlations, known as HBT interferometry, pro- 
vide a method to study the space-time structure of the emitting source [43]. 
It has been predicted that a first order phase transition would lead to unusu- 
ally large HBT-radii [44]. However, comparisons of calculations [8, 33, 45] 
with data [46, 47] have lead to the so-called HBT-puzzle: All calculations 
give a ratio of HBT-radii R ut/Rs\Ac larger than one, but the experimental 
value is of order one. The calculated values of i? Q ut are also larger and, with 
the exception of ref. [45], values of _R s ide are smaller than observed. 

It has been suggested that the solution to this puzzle lies primarily in the 
description of freeze-out [33] . This is doubtful since even calculations where 
the hadronic phase is described using a transport model cannot describe the 
data correctly [45], even though this kind of freeze-out description should 
be more reliable. Some other theoretical uncertainties, such as the order of 
the phase transition and choice of freeze-out temperature, were addressed 
in ref. [8]. As shown in fig. 7, it was found that the freeze-out temperature 
has quite a small effect on the radii nor were any of the tested EoSs (EoS 
with strong first order, weak first order and cross-over phase transitions) 
able to provide an acceptable reproduction of the data. Since the EoS with 
a smooth cross-over phase transition (CI in fig. 7) is closest to the data, it 
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Fig. 7. Rout/Rside as function of k t for Tf = 130 (left panel) and 80 MeV (right 
panel) using different equations of state compared to STAR data [46]. The calcu- 
lation and figures are from ref . [8] . 



is possible to claim that HBT measurements favour an EoS with smooth 
cross-over. This is particularly interesting since, as argued in section 3.3, 
elliptic anisotropy seems to favour a moderately strong first order phase 
transition. The explanation to this seemingly contradictory behaviour, as 
well as to the entire HBT-puzzle, is still unknown at present. 



4. Summary 

Hydro dynamical models have been very successful in explaining the sin- 
gle particle RHIC data at low pt. The pt spectra and anisotropies in central 
and semicentral collisions are well reproduced for pt < 1.5-2 GeV and the 
P/tt ratio at pt ~ 2 GeV/c has a simple explanation due to flow. Especially 
impressive has been how hydrodynamics is able to create simultaneously el- 
liptic anisotropy of negative hadrons which is large enough and anisotropy 
of protons which is small enough to fit the data. If one considers solely this 
data the collision system behaves like a thermal system. 

However, the reproduction of the HBT-radii has been unsuccessful so 
far. It is unclear whether one should refine the final freeze-out process, 
hadronization process, or initial state to reach an acceptable description of 
the data. Especially puzzling is the fact that the HBT-radii seem to favour 
a relatively stiff equation of state with a crossover phase transition, whereas 
elliptic anisotropy of protons favours a soft equation of state with a first 
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order phase transition. 
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